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We extract the leading effective range corrections to the equation of state of the unitary Fermi gas 
from ab initio Fixed-Node Quantum Monte Carlo (qmc) (fnqmc) calculations in a periodic box using 
a Density Functional Theory (dft), and show them to be universal by considering several two-body 
interactions. Furthermore, we find that the dft is consistent with the best available unbiased qmc 
calculations, analytic results, and experimental measurements of the equation of state. We also discuss 
the asymptotic effective-range corrections for trapped systems and present the first qmc results with 
the correct asymptotic scaling. 
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T-he fermion many-body problem plays a funda- 
mental role in a vast array of physical systems, from 
dilute gases of cold atoms to nuclear physics in nuclei 
and neutron stars. The universal character of this prob- 
lem - each system is governed by a similar microscopic 
theory - coupled with direct experimental access in cold 
atoms, has led to an explosion of recent interest (see 
Refs. [1] for reviews). Despite this broad applicability, 
we still do not fully understand even the simplest sys- 
tem: the "unitary gas" comprising equal numbers of 
two fermionic species with a resonant s-wave interaction 
of infinite scattering length q s — > 00. Lacking any scale 
beyond the total density n + = n Q + n b , the unitary gas 
admits no perturbative expansion and requires experi- 
mental measurement or accurate numerical simulation 
for a quantitative description. Typical Quantum Monte 
Carlo (qmc) calculations, however, can access at most a 
few hundred particles, while experiments can measure 
only a handful of properties. Density Functional The- 
ory (dft) provides a complementary approach through 
which one may extrapolate these results to large systems 
beyond the reach of direct simulation. The question 
of how the unitary gas approaches the thermodynamic 
limit has also been studied in [IH^J. 

In this paper, we consider the effects of a finite ef- 
fective range r e on the unitary gas. Our motivation is 
two-fold. First, neutron matter - well approximated 
by a unitary gas |6J - differs primarily due to a finite 
range. Characterizing the finite range effects therefore 
have physical relevance. Second, we wish to directly use 
a dft - a finite-range version of the Superfluid Local 
Density Approximation (slda) - to fit qmc simulations 
and extract thermodynamic properties without having 
to first extrapolate to zero range as was done in I3J. 
Directly fitting the finite range qmc data provides a 
much more stringent test of the slda. We use this finite- 
range dft to extrapolate to the thermodynamic limit the 
linear range dependence of the equation of state, and 



demonstrate its universality by simulating three differ- 
ent potentials. Two of the potentials include a repulsive 
core to address issues of contamination by deep bound 
states. We also show that the slda consistently fits all 
available unbiased zero-range ab initio results for the 
symmetric unitary gas. Finally, we present qmc results 
for trapped systems that demonstrate the correct asymp- 
totic scaling as predicted by the low energy effective 
theory for the unitary gas. 

Here we consider symmetric T = systems compris- 
ing equal numbers of two neutral Fermi species with 
equal mass with a short-range interaction. These sys- 
tems are directly realized by two of the lowest lying 
hyperfine states of 6 Li in cold atomic systems, and ap- 
proximately realized in dilute neutron-rich matter in 
the crusts of neutron stars. At sufficient dilution, the 
interaction can be characterized by the two-body s-wave 
phase shifts through the effective range expansion (see 
for example [7]) 

kcot6 k = — + I^ + (k 4 ) (1) 

where q s is the s-wave scattering length and r e is the 
effective range. The zero-range unitary limit is realized 
when the scattering length is tuned q s — > 00 and the 
system is diluted such that kr e — > 0: this is referred to 
as the symmetric Unitary Fermi Gas (ufg). 

The lack of scales implies that the symmetric ufg 
is fully characterized by the universal Bertsch parame- 
ter [8 J £,s = £/£fg/ where £fg = (3/5)n + Ep is the energy 
density of a free Fermi gas with the same total density 
n + = kp/(37t 2 ), and Ep = rL 2 k 2 /2m is the Fermi energy. 

In 6 Li cold-atom experiments (see [9J for details), 
q s « 00 can be tuned using the wide magnetic Feshbach 
resonance at 834.1 (15) G [10] with an effective range of 
r e ~ 4.7 nm, while the gas can be cooled at densities of 
1 /k F w 400 nm so that kpr e ~ 0.01 . In dilute neutron mat- 
ter a nrL w -18.9(4) fm [11] and r nTT « 2.75(1 1)fm mU, 



2 



while densities are on the order of 1 /kf ~ 1 fm: thus, 
IcpTe w 3 is several orders of magnitude larger than in 
cold-atom systems. 

Although there are formal ways of dealing with the 
divergences introduced by the zero-range limit (see [ [13) 
for an interesting approach), most ab initio calculational 
techniques require an explicit regulator in the form of 
a finite-range potential or a lattice cutoff. To extract 
the unitary parameters thus requires an extrapolation 
to zero effective range. Range effects in the ufg are 
also discussed in P4] (large r e ), [4] (Fixed-Node qmc 
(fnqmc)), and in [[15] [16] (Bogoliubov-de Gennes (Bdc) 
approximation) . 



By varying the Ansatz, we obtain an upper bound on 
the ground-state energy. 

We use the trial function introduced in [17]: 



¥ T =.A[*(r 11 04>(r22<)---4>(r r 



where A antisymmetrizes over particles of the same spin 
(either primed or unprimed) and f (r) is a nodeless Jas- 
trow function introduced to reduce the statistical error. 
The antisymmetrized product of s-wave pairing func- 
tions 4H r ij') defines the nodal structure: 

4>(r)=^a|| n ||e ik »- r +p(r). 



I. SUMMARY 

Here we present a summary of our results. We use 
a variational fnqmc algorithm to find upper bounds 
on the energy for systems of 4 to 66 particles in a peri- 
odic box for a variety of effective ranges kpr e 0.3 and 
for different potentials with the same scattering length 
q s = 00 and range. We fit these directly with a modi- 
fied dft that models the range dependence in order to 
extrapolate to the thermodynamic an upper bound on 
the Bertsch parameter £,, and the leading order universal 
effective range dependence C e : 



^(k F r e ) < C + C*l<Fre + 0(k F r e ) 2 
£,*= 0.3897(4) Ce= 0.127(4). 



(2a) 



By comparing several potentials, we confirm that these 
are indeed universal. The fnqmc results contain a sys- 
tematic error due to the variational nature of the method. 
To better understand this, we also fit with the dft a col- 
lection of unbiased exact, qmc, and experimental results 
for systems with 2 to 1 6 particles, obtaining a best fit of 



£ s = 0.3742(5). 



(2b) 



We also demonstrate for the first time, qmc results for 
trapped systems that exhibit the correct asymptotic be- 
havior in the thermodynamic limit. 



The sum is truncated (we include ten coefficients) and 
the omitted short-range tail is modelled by the phe- 
nomenological function (3 (r) chosen to ensure smooth 
behavior near zero separation. We use the same form 
for (3(r) as in IfiHl and vary the 10 coefficients oi» n » for 
each N + and for each different two-body potential to 
minimize the energy as described in Ref. [19]- The same 
ansatz suffices for different effective ranges, but an inde- 
pendent optimization is required for each N+. 
We compare the following potentials: 



V PT (r) =4fi 2 sech 2 ( K ir), 

V 2G (r) = 3.144li 2 ( e^ 2 / 4 - 4 e -^ rl ) , 



V 2 eM = 4.764li 2 e 



2e 



-2\xr 



(4a) 

(4b) 
(4c) 



These potentials are all tuned to have infinite two-body 
s-wave scattering length. The first potential (4a! is of 
the modified-Poschl-Teller type; the second \^b\ and 
third fcc\ potentials have a repulsive core. When tuned 
to unitarity, the effective range r e is proportional to 
as shown in figure [1] 

One criticism of purely attractive potentials - includ- 
ing the widely used modified-Poschl-Teller potential fea\ 
- is that they may contain deeply bound states where 
many particles lie within the range of the potential. For- 
mally, the ground state is thus not the universal dilute 



II. QMC MODEL 

We use a fnqmc algorithm to simulate the Hamiltonian 

where V(r) is an inter-species interaction (off-resonance 
intra-species interactions are neglected). The fnqmc 
algorithm projects out the state of lowest energy from the 
space of all wave functions with fixed nodal structure as 
defined by an initial many-body wave function (ansatz). 
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Figure 1. (color online) Finite range potentials Q used in the 
Hamiltonian {3} for our qmc bounds. 
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Figure 2. (color online) Effective range dependence of the 
ground-state energy-density £,(k F r e ) = £/£ F g of N + = 66 
fermions in a periodic cubic box in the unitary limit. The 
points with error-bars are the raw qmc results and the bands 
are the 1 a error bands of polynomial fits. The upper (blue) 
curve that extrapolates to £. = 0.3898(5) is the new quadratic 
fit to the modified-Poschl-Teller potential {4a} . The middle 
(green) curve that extrapolates to £, = 0.3885(5) is the quadratic 
fit for the new double-Gaussian potential {4b). Finally the 
lower (cyan) curve that extrapolates to £, = 0.3902(7) is the 
quadratic fit to the double-exponential potential I pjr) . 

ufg, but some tightly bound state that is highly sensitive 
to the range. In principle, this state may contaminate 
the variational qmc calculation, but in practice, there is 
insufficient overlap between the variational wave func- 
tion and this deep bound state. (Simulations longer by 
several orders of magnitude would be required to see 
the influence of such low-energy states.) 

The repulsive cores of I4E} and ( [4c) help allay these 
concerns by reducing the possibility of contamination 
from deeply bound states. We find agreement between 
the purely attractive Poschl-Teller potential and these 
repulsive potentials, demonstrating that all three poten- 
tials may be used to calculate properties of the ufg, and 
verifying the model-independence of the universal pa- 
rameters. We show the upper bounds for the energy 
of N+ = 66 particles at various effective ranges in fig- 
ure E] For ranges less than k F r e < 0.3 a three-parameter 
quadratic model is sufficient to extrapolate to zero range 
without a systematic bias. This fit is shown in table [I] for 
the three potentials, and the magnitude of the quadratic 
parameter can be used to estimate the linear regime. 

In [3], each N + was independently extrapolated to 
zero effective range, then the unitary slda dft was fit 
to the extrapolated results. It was claimed that a cubic 
fit was required to extrapolate the results for k F r e < 0.35 
to zero range, however, the smallest ranges k F r e < 0.1 
had a small systematic bias in the energy due to the 
Trotter decomposition of the many-body propagator 
e -HST _ e -VST/2 e -KST e -VST/2 + 0( v 5t )3 w h e re 5T/h 
is the imaginary time-step. Since the potentials Q scale 



roughly as V oc u 2 cx r^ 2 , for small ranges, one needs a 
very small imaginary time-step, which is computation- 
ally expensive. The extrapolated values of £, were only 
underestimated for the larger systems (by ~ 3%), but 
extracting the slope of £,(k F r e ) requires higher accuracy. 
Here we have carefully simulated with smaller time- 
steps (for the ranges considered here, 5tF_p w 5 x 10~ 6 is 
sufficient to avoid any bias) to find that, for kpr e ;$ 0.3, 
a quadratic (but not linear) fit is sufficient. We also no 
longer use an independent zero-range extrapolation for 
each N + . Instead, we use a generalized finite-range- 
Slda to fit all of the finite-range-QMC results with a 
common set of parameters. This requires simultane- 
ous consistency over all ranges and all particle numbers, 
providing a more rigorous test than independently ex- 
trapolating each N + . 

III. SLDA DFT WITH FINITE RANGE 

As was shown in [j3j, the finite-size ("shell") effects 
in f,s(N + ) can be well modelled by a simple local dft 
for the unitary Fermi gas, but are not even qualitatively 
reproduced by adding only gradient or kinetic correc- 
tions l2iT|23~| . In this paper we retain the same three- 
parameter form originally introduced in Ref. [241 (called 
the slda), but present a simple generalization that ac- 
counts for finite-range effects. With this generalized 
form, we can directly fit the qmc results without the 
need to first extrapolate to zero range. We first briefly 
review the form of the slda dft, then discuss the finite- 
range generalization. 

The slda dft is formulated in terms of three local 
densities (see [25] for a review): the total density n + , the 





£.66 (= a ) 


Ce(66) (= ttl ) 


a 2 


x r 2 


V PT 


0.3898(4) 


0.14(1) 


-0.07 


0.20 


v 2G 


0.3885(4) 


0.14(1) 


-0.08 


0.40 


V 2E 


0.3902(5) 


0.12(1) 


-0.03 


0.30 



Table I. Comparison of the zero-range extrapolations of 
£(k F r e )/£ FG = £, 66 + C e (66)kpr e + Q 2 (k F r e ) 2 + 0(r|) for 
quadratic fits of the N+ = 66 qmc results. These values are 
higher than, but consistent with the value t e (66) = 0.11(3) 
reported in [20). The extrapolations of these parameters to the 
thermodynamic limit N+ = 00 are listed as £, = a and £ e = 
in the £, block of table [n] We include the quadratic coefficient 
simply to show that the qmc can be fit using a linear form for 
k F T e < e a t, s |ai/a2l to an absolute accuracy of about e aF , s : we 
do not have any a priori reason to believe that this parameter 
is universal. The systematic error due to neglecting the cubic 
terms is on the same order as the quoted 1 cr statistical errors. 



4 



total kinetic density x+, and an anomalous v: 

n+ =2^|v n | 2 ~(ata) + (btb), 
n 

t+ = 2^JVv n | 2 ~ (Vat ■ Va) + (V6 f ■ Vb), 

n 

v = ^u n v^ ~ (ab). 
n 

These are expressed in terms of the Bogoliubov quasipar- 
ticle wave functions u n (r) and v n (r) - sometimes called 
"coherence factors". 

The three-parameter slda may then be expressed as 



£ SLDA = ^ (f T + + p A( 37 r 2 ) 2 / 3 n 5 + /3 J + gvtv, 

where oc = m/m e ff parametrizes the inverse effective 
mass; |3 parametrizes the self-energy; and g parametrizes 
the pairing interaction. In the presence of pairing, the 
local kinetic and anomalous densities are divergent 

-> A 

lim v(x,x + 5) -> — - + v r (x) + 0(5), 
5^0 5 



lim t + (x, x + 5 ) 

6^0 



+ T r (x) + 0(5), 



where A v , A T , v r and x r are finite. One must regulate 
the theory if one wishes to maintain a local formulation, 
which greatly simplifies the computational aspects of 
the dft. The most general form of a local functional 
involving these three densities is a function of these four 
finite quantities, but restricting the form to bounded 
functionals is somewhat non-trivial [26], and we shall 
not consider these generalizations here. 

We note that the 1 /5 divergence corresponds to a long 
1/k 2 momentum tail in the Fourier transform of the 
anomalous and kinetic densities. This follows from the 
short-range nature of the potential as has been empha- 
sized by Tan [13]. The most straightforward route is to 
simply introduce a momentum cutoff k < k c and then 
define the theory in the limit of large cutoff. The local 
densities then behave as 



t + = A T A + T r + 0(A 



A v A + v r + 0(A- 



where A = j k~ 2 d 3 k/(27t) 3 = k c /27t 2 . Within the single- 
particle framework of the dft, these are related to the 
gap A: A x = 2m|A| 2 /a 2 , and A v = A/a. Similar short- 
range behavior is expected in the physical density distri- 
butions where the coefficients A x and A v are related to 
the Tan's "contact" C - for example, A v = v / 2C/2m - and 
it is tempting to interpret \/2C « 2mA/a as a prediction 
of the dft, especially at unitarity where they seem to be 
related numerically. This cannot hold in general: in par- 
ticular, the contact C is related to the short-range nature 
of the interaction and persists in the normal phase (ei- 
ther meta-stable or above the critical temperature T > T c ) 



where the order parameter A vanishes [271. The 
coupling constant may be expressed 



inverse 



,-1 



:n! /3 /T-A/a, 



where y is the third dimensionless parameter character- 
izing the slda. 

The equations of motion follow by minimizing the 
total energy E = Jd 3 x £slda with respect to the oc- 
cupation factors u and v subject to the constraints of 
fixed total particle number N + and normalization. This 
leads to the following single-particle Hamiltonian for 
the Bogoliubov quasiparticle wavef unctions: 



K At 
A -K 



-VaV 
2m 



\i + U 



where U = 3£/3n + , and A = — gv. These must be solved 
self-consistently to find the stationary configurations. 
With infinite cutoff, the self-consistency equations be- 
come 



U = |3F_ F 



|A|< 



1 3/2 ' 

3n / y 



n 



1/3' 



The mean-field BdG equations may be recovered by set- 
ting a = 1, |3 = 0, and replacing n!/ 3 /y = An/ a. The 
resulting functional contains no explicit density depen- 
dence, and so contains no self-energy U = 0. The slda 
differs from the BdG equations by the inclusion of an 
effective mass and a self-energy. 

The slda functional is defined by the three dimen- 
sionless constants oc, P, and y. In practice, we use the 
homogeneous solution to the gap equation in the ther- 
modynamic limit to replace p and y with the more 
physically relevant parameters 



a : 



m 

rn-eff' 



£ 



A 



as discussed in detail in appendix |A| [see Eq. A4I. 

To extend the functional to finite range, we simply let 
the three parameters a, £,, and r| depend on the dimen- 
sionless combination kpr e . This introduces an additional 
explicit density dependence in the functional through 
kp oc n^ 3 and the self-energy must be modified ac- 
cordingly. The use of the nonlinear relationships | |A4[ 
between the polynomial form for a(kpr e ), r|(kpr e ), and 
£,(kpr e ) and the parameters of the function makes this 
complicated to write down, but numerically it is straight- 
forward to propagate these derivatives using, for exam- 
ple, automatic differentiation tools such as theano II28I . 

For the small ranges considered in this paper, we find 
that a quadratic parametrization suffices: 

a,£,,r| = qo + Qik F r e + a 2 (k F r e ) 2 . 

(Including higher order terms leads to no significant 
improvement in the quality of the fits.) This finite-range- 
Slda thus has 9 independent parameters - the three 
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coefficients a n for each of the parameters a, £,, and n. In 
comparison, the procedure of independently extrapolat- 
ing each N + to zero range introduces 3 new parameters 
for each N + in addition to the three slda parameters, 
effecting a significant increase in the total number of 
fitting parameters. Note also that the new fits directly 
use the Qmc results - including their sub-percent statis- 
tical errors - rather than the extrapolated error bar from 
zero-range extrapolation: thus the new fitting procedure 
places the slda under a significantly more stringent test. 

We only expect this extension to model the effective- 
range dependence in universal regions. In particular, 
a true finite-range interaction would naturally regulate 
the system, eschewing the need for an additional cutoff 
in the dft. For example, in the mean-field approxima- 
tion, the use of a finite-range separable potential gv^Vq 
with decaying form-factors gives rise to a momentum- 
dependent gap Ak oc Vj< regulating the anomalous den- 
sity at large momenta. Introducing such a natural reg- 
ulation into the dft, however, will likely require the 
introduction of some form of non-locality, which sig- 
nificantly complicates the computational aspects of the 
theory. 

In principle, one could also introduce a dependence 
on the scattering length a and temperature T in a similar 
manner, making the coefficients functions of kpa and 
kpjT/Ep respectively. Unlike the case with the effective 
range, however, the dependence on these parameters 
must be modelled for all values since the unitary limit 
corresponds to kpa = ±00 and T/Ep = 0, while for finite 
a and T, the zero-density limit (at the edge of a trapped 
cloud for example) is described by kpa = and T/Ep = 
00; hence, any physical system close to unitarity explores 
virtually all values of these functions, requiring a careful 
and complete characterization. 



IV. RESULTS 
A. Box 

The results of this 9-parameter fit to the qmc data-points 
with effective ranges 0.03 < kpr e ^ 0.33 are shown in 
table [n] The fit to 60 points with 4 sC N+ $: 130 for the 
V2G potential has a reduced chi squared Xr = 5. The 
fit to 70 points for 4 < N + < 130 to the V P y potential 
has Xr = 7- We suspect that this is due to approximat- 
ing the effective range dependence with a purely local 
functional as discussed earlier. 

As before [3], the best fit gap parameter ti and inverse 
effective mass oc are inconsistent with the values r\ = 
0.50(5) and a = 1.09(2) obtained from the N + = 66 qmc 
quasiparticle dispersion relation P9] [30] , and the values 
r| = 0.45(5) IfJiJ and n = 0.44(3) [[ptj extracted from 
experimental data. As we shall see below (see Eq. [7}, 
this may be due to the fixed-node approximation. 
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£,PT 


0.3911(4) 


0.111(3) 
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£.2G 


0.3900(3) 


0.111(2) 
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0.90(1) 


-0.85(7) 








87^81 


—0.82(4) 






«PT 


1.303(10) 


—0.71 (8) 






a 2G 


1 .289(7) 


—0.69(3) 






Quadratic 






a 2 


X 2 r 


£pt 


0.3903(7) 


0.121(10) 


0.00(3) 
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£.2G 


0.3890(4) 


0.128(4) 


-0.06(1) 


5 


npT 


0.99(3) 


-2.1(4) 


3(1) 




n2G 


0.879(7) 


-0.84(3) 


0.00(3) 




«PT 


1.34(2) 


-1.6(4) 


5(2) 




<*2G 


1.292(7) 


-0.73(6) 


0.1(2) 





Table II. Best fit slda parameters for linear (quadratic) 6- 
parameter (9-parameter) models: coefficients ao, a-), (and a 2 ) 
for each parameter £,, r\, and a. Note that the parameters a and 
r| should be positive, requiring positive higher-order terms that 
are not properly constrained by our qmc which only simulates 
kpT e < 0.3: larger ranges require higher-order polynomials (or 
a different model). 



Werner and Castfn [33! showed that the many-body 
energy density depends linearly on the effective range 
in the zero-range limit 



£fg 



£,s + C e kpr e + 0((kpr e 



(6) 



where the coefficient £ e is a universal constant within 
Galilean invariant continuous space models. The value 
for this coefficient was first estimated £ e = 0.046(7) [34! 
by fitting |6} to the exact two-particle solution in a trap. 

The value for this coefficient for N = 66 particles 
C e (66) = 0.11(3) was calculated using Auxiliary Field 
qmc (afqmc) in JSp] (see table [I] for comparison) and is 
likely independent of other universal parameters such 
as £, or the contact C pTj . The finite-range-SLDA allows 
us to extrapolate this result to the thermodynamic limit 
(see table [II} where we find £s = 0.127(4) by averaging 
the linear £, coefficient a, for both V PT (C e = 0.121(10)) 
and V 2 g (Ce = 0.128(4)) results. Note that these are con- 
sistent, demonstrating the universality of this coefficient. 

Unfortunately, since the fnqmc can only provide an 
upper bound on the energy, £, is systematically overesti- 
mated due to the nodal constraint. An improved nodal 
structure would lower all energies, however, and there 
is no a priori reason to suspect as large a bias for Ce- 

To address the potential systematic error introduced 
by the fixed-node approximation, we apply the same 
analysis to the recent unbiased calculations and mea- 
surements shown in table III These include zero-range 
extrapolations of two exact diagonalizations for N + = 
4 [35I, zero-range extrapolations of afqmc results for 
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Figure 3. (color online) Comparison of slda fits at zero range with zero-range extrapolated qmc upper bounds (blue) with all 
unbiased zero-range extrapolations (green) from I20I [35} listed in table [Hi] The light (yellow) band is the experimental value 
of £,s [9]. In addition, we fit the exact £2 = —0.4153 • • • value discussed in appendix [B] (not shown in the plot). Note that this 
comparison allows one to assess the fnqmc bound, which is tight for N + < 6. 



N+ = 4 |35j and for N+ e {4,14,38,48,66} J20], and 
experimental measurements of 6 Li for N + w 10 6 B9J. (Al- 
though not strictly at zero-range, the error induced by 
the non-zero range in the 6 Li experiments should be less 
than 0.003 (see also section |TVB] ).) 

We use these points to fit our three-parameter zero- 
range slda, finding: 

£,5 =0.3742(5), a =1.104(8), r, = 0.651 (9). (7) 

These error estimates must be taken with a grain of salt 
since not all of the error bars quoted in table III are 
1 ex normal standard deviations. This is reflected in the 
rather small xj = 0.2 of the fit. The results of this full fit 
are shown in figure 3] 

This partially addresses the suspiciously large value 
of n found by fitting fnqmc results (see table It ap- 
pears that a large part of the previous discrepancy is 
due to the fixed-node approximation which works well 
for small systems, but systematically overestimates the 
energy of large systems. (The variational wavef unction 
has the same number of parameters for all system sizes, 
we therefore expect it to better match the simpler nodal 
structure of small systems than the more complicated 



nodal structure of larger systems.) The gap n still ap- 
pears to be too large, which may be a problem when one 
tries to fit odd systems. 



N+ 




Method 




2 


-0.415332919- •• 


exact (see section 


B 


4 


0.288(3), 0.286(3) 


exact diagonaliza 


tion|35j 


4 


0.28(1) 


afqmc I35 




4 


0.280(4) 


AFQMC |20 




14 


0.39(1) 


AFQMC Ho 




38 


0.370(5), 0.372(2), 0.380(5) 


AFQMC Ho 




48 


0.372(3), 0.367(5) 


AFQMC |20 




66 


0.374(5), 0.372(3), 0.375(5) 


AFQMC |20 




10 6 


0.376(5) 


experiment |5j 



Table III. Unbiased zero-range box energies. Most are extrap- 
olated afqmc results except as noted. The £4 values are con- 
sistent with our upper bounds 0.2839(3) (V PT ), and 0.2829(3) 
(V2g)- This agreement indicates that the systematic error due 
to the fixed-node constraint is sub-percent for N+ = 4. 

The results shown in Fig. [2] may help understand the 
finite-size effects seen in neutron matter and neutron 
drops. In neutron matter at kpa = —10 the difference 
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Figure 4. (color online) Ground-state energy of the harmoni- 
cally trapped unitary Fermi gas (in units where hcu = 1 ) scaled 
to demonstrate the asymptotic form jgf predicted by the low- 
energy effective theory of Ref. 1 3S | . The slda with quadratic 
fit in table [II] (dashed blue line) and unbiased fit {7) (dotted 
blue line) is compared with zero-range results for N + e {4, 6} 
from Ref. [39I (black xs), and finite-range qmc results from 
Ref. I40I (upper red dots) and Ref. [41] (middle green pluses). 
The latter have significantly lower energy, despite having a 
slightly larger effective range, suggesting that the wave func- 
tions in Ref. [40! were not fully optimized. A more thorough 
optimization and extrapolation to zero effective-range yields 
the lowest points (cyan dots) which exhibit the correct scaling 
at large N + , approaching the thermodynamic value of £,. We 
also include at x = the fit £, from Eq. ^ (cyan diamond) and 
the light (yellow) experimental band [9]. 



between the energies of N + = 20 and N + = 44 particles 
is roughly 12% I36J. The range of kpr e values shown in 
Fig. [2] is too limited too allow an accurate extrapolation 
to nuclear ranges. Even so, simple extrapolations of the 
energies of N + = 14 and N + = 38 particles to kpr e = 1 .45 
using linear and quadratic forms lead to shell effects on 
the order of 10-20%, which is consistent with the (finite 
scattering-length) results seen in neutron drops [37] . 



B. Harmonic Traps 

As an application, we show here how the universal 
effective range dependence ^ affects the energy of 
particles in an isotropic harmonic trapping potential 
V(r) = mu) 2 r 2 /2 using the Thomas-Fermi (tf) approxi- 
mation. The local chemical potential is \i[r) = |xq — V(r), 
and the equation of state 



£,S + ^CekFTe + ■ 



thereby establishes the local density and energy-density 
within the tf approximation out to the maximum tf 



radius of R = ^/ly^/m/u). Including these first two 
terms we thus obtain 



N + 
E 

hcu 



u) 3 R 6 
24^372 

oj 4 R 8 
64f, 3 / 2 



32cu 4 R 7 
; 175tt£, 3 
16cu 5 R 9 
225tt£, 3 



+ ..., 



In the zero-range limit, the energy of a trapped unitary 
gas may be calculated using the low-energy effective 
theory I38J and has the form 



hcu-(3N_ 
4 



^4/3 



■6\Z27T 2 £,(2 Cl -9c 2 )(3N + )- 2/3 + 0(N 



-7/9 



where the leading order term is the well-known tf ex- 
pression (see for example [42]). The next-to-leading 
order term is directly related to the q 2 coefficient of the 
static-response and the coefficients have been estimated 
using the e-expansion [22]. The asymptotic corrections 
are due to boundary effects beyond the validity of the 
effective theory. 

This naturally suggests the introduction of the param- 
eter x = (3N + )~ 2 / 3 so that the asymptotic behavior of E 
is linear in x. The square of the energy E 2 also exhibits 
linear asymptotic behavior, 



16E 2 



h 2 cu 2 (3N+) 8 / 3 



£, + cx + 0(x 



7/6i 



(9) 



where c = — 12v / 27t 2 £, 3 / 2 (2ci — 9 02)- We prefer this form 
as f, appears as the intercept and note that the relation- 
ship is remarkably linear, as can be seen in figure [4] 

It is interesting that, in the non-interacting system, 
shell-effects appear at the same linear order x, leading 
to a fundamental uncertainty in the coefficient 0.67 < 
c < 1.7. Pairing suppresses these shell effects, and they 
are virtually non-existent in the unitary gas leading to a 
well-defined value of c. Note that the tf approximation 
contains only the leading order term: i.e. c = 0. 

In the tf approximation, the leading-order effective- 
range correction Cekpr e leads to a super-leading order 
(in N + ) correction to ((9): 



16E' 



h 2 cu 2 (3N + ) 8 



/3 



^+1.17^ 7 fv / tux- 1 /4 + (r 2 ). (10) 



(The coefficient is 1.17 = 2 25 / 2 /1575tt.) The singular 
x~ 1 / 4 demonstrates that, as N + gets large, the central 
density becomes large and kpr e corrections play an in- 
creasingly significant role. The analysis is therefore only 
valid in a limited regime where the system is sufficiently 
large that the tf approximation is valid, but where the 
central density is small enough that kpr e remains small. 
It illustrates how a finite effective range will alter the 
linear asymptotic behavior expected in figure [4] 
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In figure [4] we show new fixed-node qmc results that 
have been extrapolated to zero-range using a quadratic 
polynomial in kpr e . These results represent the first 
ab initio calculations to demonstrate the correct linear 
asymptotic scaling as predicted by the effective theory. 
In particular, all previous results start to "turn up" as 
they approach the thermodynamic limit. While this 
is qualitatively consistent with the expected divergent 
x _1 ' 4 behavior expected of a finite-range, the effect does 
not agree quantitatively: eq. [10] predicts the divergence 
to set in at a larger N + than seen in figure [4] We suspect 
that the incorrect scaling of previous trapped results 
indicates the presence of spurious length scales (but 
in principle, could also signify spurious breaking of 
another symmetry). 

Allowing a more flexible variational wavefunction 
(green pluses I41I) improves the bound compared with 
the red dots of [40]. This seems sufficient for small sys- 
tems as witnessed by the agreement with the N + € {4, 6} 
results of [39I, but does not provide the correct asymp- 
totic behavior in larger traps where the density and pair- 
ing correlations differ substantially between the center 
and edges of the trap. To obtain the correct asymptotic 
behavior here, we include an explicit dependence on the 
center-of-mass coordinate of each pair in the variational 
pairing wavefunction (cyan dots). 

The linear scaling of our new results indicates that this 
nodal approximation does not introduce any spurious 
length scales, however, even with this extra freedom, the 
variational bound provided for trapped systems is not 
as tight as it is for homogeneous matter, and the cyan 
dots extrapolate to a somewhat higher bound for the 
value of £, « 0.4. As with the homogeneous systems, we 
find the same trend that the variational bound is tight 
for small systems, but is less accurate for larger systems 
where pairing correlations become more significant. 

Finally, we have included the slda predictions in the 
figure E] (blue curves) but do not use the slda to fit 
the results since we have not included any gradient cor- 
rections. By construction, the slda extrapolates to the 
thermodynamic value of £, used in the parametrization, 
but the slope is sensitive to the leading order gradient 
corrections that we have neglected in this paper since 
they do not contribute to homogeneous matter. This 
plot contains within it hints as to the nature of the gradi- 
ent corrections to the slda, but quantitative statements 
require further analysis beyond the scope of this paper. 



V. SUMMARY AND CONCLUSIONS 

In this work we have extensively analyzed the ground- 
state energy of strongly interacting atoms for finite effec- 
tive ranges. We present new Fixed-Node qmc results for 
inter-atomic potentials that also contain repulsive cores: 
these new potentials yield results that are statistically 
consistent with the purely attractive (modified Poschl- 
Teller) potential used in earlier works, demonstrating the 
universality of the leading finite-effective-range depen- 
dence, and addressing concerns about contamination of 
the fnqmc energies by deeply bound many-body states. 

To model these results in a common framework, we 
have minimally extended the Superfluid Local Density 
Approximation Density Functional Theory to directly 
fit the finite-range fnqmc results. Although this simple 
generalization of the slda is not completely consistent 
with the fnqmc results, it still proves to be a useful 
tool for extrapolating finite-size results to the thermody- 
namic limit. To assess the accuracy of the variational up- 
per bound provided by the fnqmc results, we have also 
fit the slda to unbiased (non-variational) exact, qmc, 
and experimental results from the literature to produce 
a working slda for modeling physical systems. This fit 
demonstrates that the three-parameter zero-range slda 
is consistent with the unbiased results. 

Finally, we have presented new qmc results for zero- 
range trapped systems. These results demonstrate, for 
the first time, the correct asymptotic behavior in the 
thermodynamic limit as predicted by the low-energy 
effective theory. 
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Appendix A: Homogeneous solutions of the Superfluid 
Local Density Approximation 

In this appendix, we describe some properties of homo- 
geneous solutions to the slda functional, both in the 
periodic box, and in the thermodynamic limit of infinite 
matter. (When these equations are applied locally at 
each point in a slowly varying external potential, one ob- 
tains the tf approximation.) As discussed in the text, we 
use the thermodynamic solutions to express the param- 
eters p and y in terms of the more physically relevant 
quantities £, and r). 

We start by rotating away the phase, taking A = |A| to 
be real. We also note that the self -energy U plays no role 
in the solution of the homogeneous equations: all effects 
are absorbed into the effective chemical potential (i^g. 
One only needs to compute the self energy U to relate 
the effective chemical to the thermodynamic chemical 
potential. Thus, the homogeneous Hamiltonian is com- 
pletely parametrized by a, |j. e ff, and A. In momentum 
space, the Hamiltonian is easily diagonalized, 



cutoff on both the discrete and continuous momenta. 
The partial sums as a function of cutoff will fluctuate as 
various lattice points enter the sphere, but the magnitude 
of the fluctuations will reduce and the resulting limit 
converges. Numerically, it is favorable to sum over 
cubic shells so that the sequence of partial sums behaves 
smoothly, allowing one to accelerate the convergence. 
However, the location of the cutoff between shells must 
be fine tuned to reproduce the correct result because 
- unlike the spherical case - the fluctuations never die 
away with a cubic cutoff. 

From the integrals one can see that the effective mass 
can be scaled out to define the following finite functions: 



a' a 
A HefA 



C 

\ a a / 



a. a. J 



u. eff ,A), 



-v r (a, ^ eff ,A), 



A 2 



T r (a, ^ eff , A). 



H 



A 



A ) = ( Uk 
ah 2 k 2 



£k 



u-k 



2m 



Vk 

^eff/ 



u k 



Vk 



e£+A 2 , 



1 + 



Vk 



e + 
E + 



In this diagonal form, the density matrix p = f p (H) 
can be computed in a straightforward manner from the 
Fermi distribution function fp(H) = 1/[1 +exp(— |3H)]. 
For reference, the zero-temperature results are: 



n+(a, Heff, A) 
2m 



d 3 k 



t-(2r 
J- d 3 k k 2 
4-(2ti) 3 2m 

J- d 3 k A 
4-j2n)3 2Ek~ 



£k 
Hk 



£k 
Ek 



The notation rf^T5 represents either a discrete summa- 
tion over box momenta ki = 2mn.i/Li or the continuous 
integral Jd 3 k /(27t) 3 in the thermodynamic limit. The 
regulated quantities r r and y r follow from these by sub- 
tracting the power-law divergences (this is equivalent to 
using dimensional regularization (43]) : 



2m 



f d 3 k 
4-(2tt) 3 
f d 3 k 
4-(2ti) 3 



k 2 ^ 
2m 
A 
2Ek 



£k 

d 3 k 
(27tp 



d 3 k mA 2 



(2tt) 3 ct 2 k 2 ' 



mA 
ok 2 ' 



Note that the subtraction integrals are continuous. In 
order to implement this regularization scheme in the 
periodic box, one must use a simultaneous spherical 



One can thus deduce that, if the volume V = L x Ly L z and 
shape of the box L are held fixed, then the tf equations 
exhibit an additional invariance under scaling a, |x e ff, 
and A by the same factor 

da _ d|j. e ff _ dA 
a ^eff A ' 

Note that this does not follow from dimensional analysis 
(a is already dimensionless) and expresses a non-trivial 
property of the tf equations. These scaling relation- 
ships allow us to express everything in terms of two 
dimensionless parameters - K = r\/<x, and the total (di- 
mensionless) particle number N + = Vn+ - through the 
dimensionless functions c[K, N + ) and d(K, N + ): 



K EE 



n 



a 

c(K,N+) = C 
3K 



N+ ee n+V, 



c(K,N + )n+, 



d(K,N_ 



-FG 



K 2 d(K,N + )£ FG . 



In the T = thermodynamic limit L — > oo (N + — > oo), 
the integrals can be performed analytically (see |j43j)- 
We start by defining: 



ko 



2m|n. eff | 
ah 2 ' 



yo 



^eff 
A 2 + - 2 



^eff 



We may then express our previous results as 



3tt 2 



mko 

47th 2 



-cc 2 h 2 k 5 



H c , D 



Oi 



8m7t 2 A 2 
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where the functions h. n , h c , and h d depend only on yo, 

, , > 3yofi/2(yo) -f3/2(yo) , , , fi/2(yo) 



h d (yo) 



4 lyol 3 / 2 

f5/2(yo) -yof3/2(yo) 



lyo 



5/2 



fa(yo) 



-TtPot(-yo) 

sin(Tra) 



through the Legendre function P a (x) which satisfies 


P«fxl 



(1-x 2 )P^-2xP^ + «(a+1)P a/ 
1 



2m 



ooi 



\J 1 - 2xcu + w 2 dco. 



Noting that n + = kp/37T 2 we can identify kp = k^h^ and 

E F = h 2 k 2 /2m = rt£ /3 lc§/2m = h 2/3 |u- eff |/a. We can then 
relate X directly to yo through the monotonic function: 



K(yc 



yo 2 - 1 



«Ef |n eff |h 2/3 (y ) H 2/3 (y 



This function has the limiting behavior: 

1/6 



X 



(t^o) ' where y ~-l, 
V^H -yo) where y ~ h 



and an application of five steps of Newton's method 
using this as a guess (splitting the input at the point 
X w 1.211 292490 where these asymptotic forms meet) 
solves the inverse problem yo(X) to machine precision. 
With this conversion we can directly express 



c(K) = c N + =0O (X) 
d(K) = d N + =0O (K) 



5tt h c (y ) 
8 J/3, / 

hn. (yo) 
-5 h d (y ) 

4 ^H 5 n /3 (y )' 



allowing the parameters a, (3, and y to be computed 
from the thermodynamic values of a, £,, and T|: 



Y 
P 



5a(37t 2 ) 2 / 3 



6c{t)/<x) 
d(Ti/a)Ti 2 6ri 2 y 
a 5(3tt 2 ) 2 / 3 ' 



(A4a) 
(A 4 b) 



We use these equations to express all of our results in 
terms of the thermodynamic values of a, £,, and r\, even 
though the functional is expressed in terms of fixed 
parameters a, (3, and y. 



Appendix B: Particles in a Box 

Here we present some details about computing the en- 
ergies E of N + = N a + Nt, particles in a cubic box of 
size L 3 . There are two conventions for expressing the 



energy of a box. We use £,n + = £(N + )/£ FG where 
£(N+) = E(N + )/L 3 . All values of £, reported in this 
paper have been converted to this normalization. The 
other convention £, box = E(N + )/Efg(N + ) normalizes the 
energy with respect to the energy of N + non-interacting 
fermions in the same box (see [3] for conversion factors). 

To further constrain our fits, we include the results for 
N+ = 2. By solving the Schrodinger equation for two 
particles in a periodic box of size L 3 with the short-range 
boundary condition 

lim¥(x,x + f) oc - +kcot5 k + 0(r), 
r^o r 



one obtains 

kcot5 k 



7tL 



S(tO = lrm Y_ 



1 



nr-11 



■AnA. 



where E = k 2 /2m r is the energy in the center-of -mass- 
frame and m r = m/2 is the reduced mass of the system 
(see for example [44] and references therein). Note that 
for non-interacting particles, Epc(N + ) = 0, thus for all 
attractive interactions, E oc k 2 < 0. This poses no prob- 
lems since only k 2 enters the formulation: for example, 
kcot 5k is a series in k 2 Q. 

As before, the summation may be performed with par- 
tial sums over cubic shells: these behave smoothly and 
are amenable to series acceleration techniques (see I45J 
for example) such as the Levin transformation. 

The energies £,2 (kpr) are shown in figure [5] for the 
potentials Q. Over the ranges considered, the results 
are virtually identical. Finally, we note that the N + = 2 
solution to the slda has only the normal solution A = 0. 
Both particles enter the k = ground state which has 
zero energy, hence we can identify p(kpr) = £,2 (kpr). 
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Figure 5. (color online) Exact ground-state energy £,2(' c F r e 
the N+ = 2 system in a box for each of the potentials (4}. 
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